
%config InlineBackend.figure_format = 'retina'
%matplotlib inline    
import seaborn as sns 
sns.set(font= "Kaiti",style="ticks",font_scale=1.4)
import matplotlib
matplotlib.rcParams['axes.unicode_minus']=False # 解决坐标轴的负号显示问题
## 导入需要的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import plotly.express as px
from fbprophet import Prophet
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.metrics import *
from sklearn.svm import SVC
from sklearn.manifold import TSNE

import pyod.models as pym
from pyod.models.cof import COF
from pyod.models.pca import PCA
from pyod.models.sod import SOD
from pyod.models.iforest import IForest
from pyod.models.xgbod import XGBOD


usedatadf = pd.read_excel("附件1：赛题A数据.xlsx")
usedatadf.head()


## 数据提取,获取区内的平均用户数数据，指标1
usenumdf = usedatadf.iloc[:,0:4]
usenumdf["number"] = usedatadf["小区内的平均用户数"]
usenumdf["时间"] = pd.to_datetime(usenumdf["时间"],format = "%Y-%m-%d %H:%M")
usenumdf.head()


## 根据模型预测值的置信区间"yhat_lower"和"yhat_upper"判断样本是否为异常值
def outlier_detection(forecast):
    index = np.where((forecast["y"] <= forecast["yhat_lower"])|
                     (forecast["y"] >= forecast["yhat_upper"]),True,False)
    return index

## 定义一个输出是否为异常值的列表
def outlier_detection2(forecast):
    outlier = []
    for ii in forecast.index:
        inde = ((forecast.loc[ii,"y"] <= forecast.loc[ii,"yhat_lower"])|(forecast.loc[ii,"y"] >= forecast.loc[ii,"yhat_upper"]))
        outlier.append(inde)
    return outlier


modeldf = usenumdf
modeldf["outlier"] = None
modeldf.head()


## 检测出的异常值的数量，也会对后面第2个问题的分类模型有影响（或者异常值检测模型）
varname = "number"  ## 需要检测的变量
plotresult = True    ## 是否输出图像

## 获取所有小区的编号
np.random.seed(1234)  ## 设置随机数种子
xiaoqu = np.unique(modeldf["小区编号"])
for ii,xq in enumerate(xiaoqu):
    ## 获取一个小区的数据
    predf = modeldf.loc[modeldf["小区编号"] == xq,["时间",varname]]
    predf.columns = ["ds","y"]
    ## 对一个小区的数据建立一个时间序列模型
    model = Prophet(growth="linear",daily_seasonality = True,
                    weekly_seasonality=True,yearly_seasonality=False,
                    seasonality_mode = 'multiplicative', 
                    interval_width = 0.95,   ## 用于控制异常值的数量（0.95，0.99，0.9等）
                    )
    model = model.fit(predf)     # 使用数据拟合模型
    forecast = model.predict(predf)  # 使用模型对数据进行预测
    forecast["y"] = predf["y"].reset_index(drop = True)
    forecast["yhat"] = np.where(forecast["yhat"]<0,0,forecast["yhat"])  # 有些0数值不能粗暴的认为其就是异常值
    forecast["yhat_lower"] = np.where(forecast["yhat_lower"]<-0.5,-0.5,forecast["yhat_lower"])
    ## 查看异常值预测情况
    outlier_index = outlier_detection(forecast)
    outlier_df = predf[outlier_index]
    print("小区:",xq,"  异常值的数量为:",np.sum(outlier_index))
    ## 获取异常值列表
    outlier_index2 = outlier_detection2(forecast)
    predf["outlier"] = outlier_index2 
    ## 将是否为异常值填写到原始数据中
    modeldf.loc[predf.index,"outlier"] = predf["outlier"]
    ## 可视化预测结果
    if plotresult:
        ## 可视化异常值的结果
        fig, ax = plt.subplots()
        ## 可视化预测值
        forecast.plot(x = "ds",y = "yhat",style = "b-",figsize=(16,7),
                      label = "预测值",ax=ax)
        ## 可视化出置信区间
        ax.fill_between(forecast["ds"].values, forecast["yhat_lower"], 
                        forecast["yhat_upper"],color='b',alpha=.2,
                        label = "95%置信区间")
        forecast.plot(kind = "scatter",x = "ds",y = "y",c = "k",
                      s = 10,label = "原始数据",ax = ax)
        ## 可视化出异常值的点
        outlier_df.plot(x = "ds",y = "y",style = "rs",ax = ax,
                        label = "异常值")
        plt.legend(loc = 2)
        plt.grid()
        plt.title("小区:"+str(xq)+"时间序列异常值检测结果")
        plt.show()

Initial log joint probability = -13.8552
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       1231.64   0.000946278       50.5977      0.8269      0.8269      130   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     105       1231.81    0.00214368         98.29   2.455e-05       0.001      170  LS failed, Hessian reset 
     188       1231.95    2.4394e-05       51.6132   4.055e-07       0.001      314  LS failed, Hessian reset 
     199       1231.96     7.377e-07        67.578       0.357       0.357      329   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     234       1231.96   5.79168e-08       75.6851      0.3249           1      387   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019001   异常值的数量为: 36

Initial log joint probability = -10.1021
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       1665.67   1.51109e-05       64.5753      0.2198      0.2198      131   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     123       1665.68   6.77103e-06       74.5389   8.208e-08       0.001      204  LS failed, Hessian reset 
     144       1665.68   1.49099e-08       75.0703      0.0856       0.382      229   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019002   异常值的数量为: 38

Initial log joint probability = -16.1202
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      59       1235.61    0.00139449       154.099   1.521e-05       0.001      102  LS failed, Hessian reset 
      99       1235.91    0.00273097        70.317           1           1      161   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     111       1236.23   0.000497332       61.7732   7.306e-06       0.001      224  LS failed, Hessian reset 
     199       1236.55   1.05246e-05       70.2324      0.4468      0.4468      340   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     227       1236.55   1.01248e-06       72.1325   1.536e-08       0.001      418  LS failed, Hessian reset 
     236       1236.55   9.25902e-08        62.382      0.1723           1      430   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019003   异常值的数量为: 31

Initial log joint probability = -16.714
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       1243.84   6.84863e-06       61.7659           1           1      145   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     146       1244.18    0.00114975       94.4238     1.5e-05       0.001      271  LS failed, Hessian reset 
     199       1244.41   5.24026e-05       73.2198      0.3835      0.3835      338   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     299       1244.43    9.6596e-06       86.3621           1           1      465   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     365       1244.46   7.81983e-09       63.6025     0.06355      0.9397      553   
Optimization terminated normally: 
  Convergence detected: absolute parameter change was below tolerance
小区: 26019004   异常值的数量为: 32

Initial log joint probability = -32.8427
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       1355.85    0.00166816       83.4561           1           1      132   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     179       1356.53   0.000162473       77.5984   2.582e-06       0.001      276  LS failed, Hessian reset 
     199       1356.54   1.38495e-05       84.9198      0.6384      0.6384      301   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     218       1356.56   3.34127e-06       54.0073   4.784e-08       0.001      369  LS failed, Hessian reset 
     261       1356.57   2.10543e-05       73.9325   3.224e-07       0.001      455  LS failed, Hessian reset 
     299       1356.57   1.06392e-06       64.9042      0.7121      0.7121      506   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     311       1356.57   6.87284e-08       58.6732      0.5036      0.5036      521   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019005   异常值的数量为: 45

Initial log joint probability = -25.0319
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       1230.53   0.000591093        90.279           1           1      127   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     139       1230.72   0.000317865       52.6359   3.866e-06       0.001      269  LS failed, Hessian reset 
     199       1230.86   5.90223e-06        73.339      0.1213           1      342   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     224       1230.86   1.25289e-07       65.3868      0.2123           1      375   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019006   异常值的数量为: 46

Initial log joint probability = -9.89945
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99       1491.57   0.000907743       114.967           1           1      118   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     132       1491.85   0.000258062       75.6505   2.816e-06       0.001      194  LS failed, Hessian reset 
     199       1492.01   8.17498e-05       67.8251      0.2193      0.2193      283   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     289       1492.48    1.0099e-07       57.2788      0.3975           1      395   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019007   异常值的数量为: 33

Initial log joint probability = -10.5983
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      97       1404.78   0.000431619       89.0976   3.551e-06       0.001      205  LS failed, Hessian reset 
      99       1404.82   0.000652859       73.8344           1           1      207   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     184       1404.97   2.90576e-06       67.0254   4.203e-08       0.001      366  LS failed, Hessian reset 
     199       1404.97   1.45266e-07       62.7546      0.6752      0.6752      390   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     204       1404.97   1.79207e-08       61.0192       0.216       0.216      397   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019008   异常值的数量为: 31

Initial log joint probability = -35.0568
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      71       902.258    0.00348396       89.4195   3.588e-05       0.001      120  LS failed, Hessian reset 
      99       902.424   4.72416e-06       71.7274      0.2976      0.2976      156   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     199       902.527   1.81757e-05       74.0284      0.4554      0.4554      295   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     299       902.539   2.95059e-06       75.2342           1           1      419   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     399        902.56     2.506e-07       70.1503      0.5793      0.5793      545   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     440        902.56   3.76361e-06       82.3118   4.333e-08       0.001      641  LS failed, Hessian reset 
     461       902.561   5.02106e-08       75.0374      0.1854           1      670   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019009   异常值的数量为: 38

Initial log joint probability = -10.4101
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99        1128.6   0.000306573       83.7598      0.4665           1      139   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     199       1129.15   4.41284e-05       84.1581      0.1782           1      271   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     226       1129.15   2.55137e-08       76.2364     0.04549      0.4854      312   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
小区: 26019010   异常值的数量为: 31


## 预测结果的查看
modeldf


## 计算每个小区的异常值数量
modeldf.groupby("小区编号")["outlier"].sum().plot(kind = "bar",figsize = (17,8))

<AxesSubplot:xlabel='小区编号'>


## 针对是否时异常值区间的数据，根据看是否有几个异常值聚集到一块进行判断即可。其它的变量也可使用同样的方式进行预测与判断是否为异常值


## 针对某个小区，只使用时序数据的模型建立
cladf =  modeldf.loc[modeldf["小区编号"] == 26019005,["时间","number","outlier"]]
cladf["outlier"] = np.int32(cladf["outlier"])
cladf


## 可视化时间序列的变化情况
plt.figure(figsize=(20,10))
sns.pointplot(data = cladf,x = "时间",y = "number",hue="outlier")
plt.show()


## 数据准备，使用此否为异常值的前24个样本作为模型的特征(不同的长度会有不同的预测效果)
featurelen = 24
X = []
Y = []
for ii in np.arange(featurelen,len(cladf)):
    X.append(cladf.number.values[(ii-featurelen):ii])
    Y.append(cladf.outlier.values[ii])
X = np.array(X)
Y = np.array(Y)


## 数据降维可视化
## TSNE进行数据的降维,降维到2维空间中
tsne = TSNE(n_components = 2,perplexity =100,
            early_exaggeration =5,random_state=123) 

## 获取降维后的数据
tsne_outlier = tsne.fit_transform(X)
print(tsne_outlier.shape)
## 可视化在二维空间中的分布
plt.figure(figsize=(10,7))
sns.scatterplot(x=tsne_outlier[:,0], y=tsne_outlier[:,1],
                style = Y,hue = Y,s = 80)
plt.legend(loc = 1)
plt.grid()
plt.title("TSNE降维后数据的分布")
plt.show()

## 异常值的分布并不那么容易识别

(672, 2)


## 针对这样的数据类型，由于异常值的数量很少，所以可以看作是一个二分类问题，也可以看作是一个新的异常值检测问题


## 先使用二分类模型对其进行检测，再介绍一些异常值识别模型的使用

## 数据切分训练集和验证机
X_train, X_test, y_train, y_test = train_test_split(X, Y,test_size=0.25,random_state=2)
print("训练数据:",X_train.shape)
print("训练数据中异常值数量:",np.unique(y_train,return_counts = True))
print("测试数据:",X_test.shape)
print("测试数据中异常值数量:",np.unique(y_test,return_counts = True))

训练数据: (504, 24)
训练数据中异常值数量: (array([0, 1], dtype=int32), array([471,  33]))
测试数据: (168, 24)
测试数据中异常值数量: (array([0, 1], dtype=int32), array([156,  12]))


## 使用基于PCA方法的异常值识别
pcaod = PCA(n_components = "mle",  #自动猜测保留的主成分数量
            n_selected_components = 4, # 计算异常值得分时使用的主成分数量
            contamination = 0.1,    # 异常值所占比例
            random_state = 123)
pcaod.fit(X_train)   ## 对训练数据进行拟合
pcaod_lab = pcaod.labels_
print("在训练集上是否为异常值判断正确的精度为:\n",classification_report(y_train,pcaod_lab))
print("在训练集上是否为异常值判断正确的f1_score为:",f1_score(y_train,pcaod_lab))
## 对测试集进行预测，并计算预测的精度
pcaod_pre = pcaod.predict(X_test)
print("在测试集上是否为异常值判断正确的精度为:\n",classification_report(y_test,pcaod_pre))
print("在测试集上是否为异常值判断正确的f1_score为:",f1_score(y_test,pcaod_pre))

在训练集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.93      0.89      0.91       471
           1       0.02      0.03      0.02        33

    accuracy                           0.84       504
   macro avg       0.47      0.46      0.47       504
weighted avg       0.87      0.84      0.85       504

在训练集上是否为异常值判断正确的f1_score为: 0.02380952380952381
在测试集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.94      0.85      0.89       156
           1       0.11      0.25      0.15        12

    accuracy                           0.80       168
   macro avg       0.52      0.55      0.52       168
weighted avg       0.88      0.80      0.84       168

在测试集上是否为异常值判断正确的f1_score为: 0.15384615384615383


## SOD是一种子空间异常值检测方法
sod = SOD(n_neighbors = 20,  ## 使用K近邻查询是的近邻数量
          ref_set = 10, ## 创建参考集的共享最近邻居的数量
          alpha = 0.85, ## 选择指定子空间的下限
          contamination  = 0.1)
sod.fit(X_train)
sod_lab = sod.labels_
print("在训练集上是否为异常值判断正确的精度为:\n",classification_report(y_train,sod_lab))
print("在训练集上是否为异常值判断正确的f1_score为:",f1_score(y_train,sod_lab))
## 对测试集进行预测,并计算预测精度
sod_pre = sod.predict(X_test)
print("在测试集上是否为异常值判断正确的精度为:\n",classification_report(y_test,sod_pre))
print("在测试集上是否为异常值判断正确的f1_score为:",f1_score(y_test,sod_pre))

OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

在训练集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.93      0.90      0.91       471
           1       0.04      0.06      0.05        33

    accuracy                           0.84       504
   macro avg       0.49      0.48      0.48       504
weighted avg       0.87      0.84      0.86       504

在训练集上是否为异常值判断正确的f1_score为: 0.04761904761904762
在测试集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.94      0.71      0.81       156
           1       0.10      0.42      0.16        12

    accuracy                           0.69       168
   macro avg       0.52      0.56      0.49       168
weighted avg       0.88      0.69      0.76       168

在测试集上是否为异常值判断正确的f1_score为: 0.16129032258064516


xgod = XGBOD(n_estimators =100, ## 树的数量
             max_depth = 10,    # 最大深度
             learning_rate = 0.1, # 学习率 
             reg_alpha = 0.01, # L1范数约束的权重
             random_state = 12)
## 注意在拟合该模型时，使用了是否为异常值的标签进了监督
xgod.fit(X_train,y_train)  
xgod_lab = xgod.labels_
print("在训练集上是否为异常值判断正确的精度为:\n",classification_report(y_train,xgod_lab))
print("在训练集上是否为异常值判断正确的f1_score为:",f1_score(y_train,xgod_lab))
## 对测试集进行预测,并计算预测精度
xgod_pre = xgod.predict(X_test)
print("在测试集上是否为异常值判断正确的精度为:\n",classification_report(y_test,xgod_pre))
print("在测试集上是否为异常值判断正确的f1_score为:",f1_score(y_test,xgod_pre))

/Users/daitu/opt/anaconda3/envs/DataVis/lib/python3.8/site-packages/pyod/models/base.py:415: UserWarning: y should not be presented in unsupervised learning.
  warnings.warn(
/Users/daitu/opt/anaconda3/envs/DataVis/lib/python3.8/site-packages/xgboost/sklearn.py:1224: UserWarning: The use of label encoder in XGBClassifier is deprecated and will be removed in a future release. To remove this warning, do the following: 1) Pass option use_label_encoder=False when constructing XGBClassifier object; and 2) Encode your labels (y) as integers starting with 0, i.e. 0, 1, 2, ..., [num_class - 1].
  warnings.warn(label_encoder_deprecation_msg, UserWarning)

[15:42:04] WARNING: /Users/travis/build/dmlc/xgboost/src/learner.cc:576: 
Parameters: { "silent" } might not be used.

  This could be a false alarm, with some parameters getting used by language bindings but
  then being mistakenly passed down to XGBoost core, or some parameter actually being used
  but getting flagged wrongly here. Please open an issue if you find any such cases.


[15:42:04] WARNING: /Users/travis/build/dmlc/xgboost/src/learner.cc:1115: Starting in XGBoost 1.3.0, the default evaluation metric used with the objective 'binary:logistic' was changed from 'error' to 'logloss'. Explicitly set eval_metric if you'd like to restore the old behavior.
在训练集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       1.00      1.00      1.00       471
           1       1.00      1.00      1.00        33

    accuracy                           1.00       504
   macro avg       1.00      1.00      1.00       504
weighted avg       1.00      1.00      1.00       504

在训练集上是否为异常值判断正确的f1_score为: 1.0
在测试集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.94      0.99      0.97       156
           1       0.60      0.25      0.35        12

    accuracy                           0.93       168
   macro avg       0.77      0.62      0.66       168
weighted avg       0.92      0.93      0.92       168

在测试集上是否为异常值判断正确的f1_score为: 0.35294117647058826


## 介绍使用SVM分类器建立分类模型（这里只介绍这一种方法，其它的就不再介绍了）
## 建立非线性SVM模型，使用rbf核
rbfsvm = SVC(kernel  = "rbf",gamma=0.0001, ## RBF核和对应的参数
             random_state= 1,C=1000)
## 训练模型
rbfsvm.fit(X_train,y_train)
## 计算在训练集和测试集上的预测精度
rbfsvm_lab = rbfsvm.predict(X_train)
rbfsvm_pre = rbfsvm.predict(X_test)
print("训练集预测精度:",accuracy_score(y_train,rbfsvm_lab))
print("测试集预测精度:",accuracy_score(y_test,rbfsvm_pre))
print("训练集预测精度:",f1_score(y_train,rbfsvm_lab))
print("测试集预测精度:",f1_score(y_test,rbfsvm_pre))
print("在训练集上是否为异常值判断正确的精度为:\n",classification_report(y_train,rbfsvm_lab))
print("在测试集上是否为异常值判断正确的精度为:\n",classification_report(y_test,rbfsvm_pre))

训练集预测精度: 1.0
测试集预测精度: 0.9285714285714286
训练集预测精度: 1.0
测试集预测精度: 0.5
在训练集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       1.00      1.00      1.00       471
           1       1.00      1.00      1.00        33

    accuracy                           1.00       504
   macro avg       1.00      1.00      1.00       504
weighted avg       1.00      1.00      1.00       504

在测试集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.96      0.96      0.96       156
           1       0.50      0.50      0.50        12

    accuracy                           0.93       168
   macro avg       0.73      0.73      0.73       168
weighted avg       0.93      0.93      0.93       168


usedatadf.columns

## 这里我们添加几个作为示例，还可使用更多的其它特征
myaddusefeature = ['上行可用的PRB个数', '下行可用的PRB个数',
       '上行PhysicalResourceBlock被使用的平均个数', '下行PhysicalResourceBlock被使用的平均个数',
       '上行PUSCH的PhysicalResourceBlock被使用的平均个数','ERAB建立成功总次数', 'ERAB建立尝试总次数',
       'ERAB异常释放总次数', 'ERAB正常释放总次数', '系统间切换出ERAB正常释放总次数', 'eNodeB内同频切换出成功次数',
       'eNodeB间同频切换出成功次数', 'eNodeB内同频切换出执行次数', 'eNodeB间同频切换出执行次数',
       'eNodeB内异频切换出成功次数', 'eNodeB间异频切换出成功次数']


## 针对某个小区，只使用时序数据的模型建立
cladf2 =  usedatadf.loc[usedatadf["小区编号"] == 26019005,myaddusefeature]
cladf2


cladf2[(25-featurelen):25].apply(func = np.mean,axis= 0)

上行可用的PRB个数                                100.000000
下行可用的PRB个数                                100.000000
上行PhysicalResourceBlock被使用的平均个数            10.882679
下行PhysicalResourceBlock被使用的平均个数            26.663988
上行PUSCH的PhysicalResourceBlock被使用的平均个数       6.953225
ERAB建立成功总次数                              6219.958333
ERAB建立尝试总次数                              6223.375000
ERAB异常释放总次数                                 0.458333
ERAB正常释放总次数                              6287.125000
系统间切换出ERAB正常释放总次数                           0.000000
eNodeB内同频切换出成功次数                          376.875000
eNodeB间同频切换出成功次数                          844.708333
eNodeB内同频切换出执行次数                          377.000000
eNodeB间同频切换出执行次数                          845.458333
eNodeB内异频切换出成功次数                            0.208333
eNodeB间异频切换出成功次数                           49.375000
dtype: float64


## 数据准备，使用此否为异常值的前24个样本作为模型的特征(不同的长度会有不同的预测效果)
featurelen = 24
X2 = []
for ii in np.arange(featurelen,len(cladf2)):
    X2.append(cladf2[(25-featurelen):25].apply(func = np.mean,axis= 0).values)
X2 = np.array(X2)
print(X2.shape)

(672, 16)


## 将X2和前面的X拼接，获得新的X
X = np.hstack((X,X2))
print(X.shape)

(672, 40)


Y.shape

(672,)


## 对拼接后的X进行标准化处理
Xs = StandardScaler().fit_transform(X)
## 数据切分训练集和验证机
X_train, X_test, y_train, y_test = train_test_split(Xs, Y,test_size=0.25,random_state=2)
print("训练数据:",X_train.shape)
print("训练数据中异常值数量:",np.unique(y_train,return_counts = True))
print("测试数据:",X_test.shape)
print("测试数据中异常值数量:",np.unique(y_test,return_counts = True))

训练数据: (504, 40)
训练数据中异常值数量: (array([0, 1], dtype=int32), array([471,  33]))
测试数据: (168, 40)
测试数据中异常值数量: (array([0, 1], dtype=int32), array([156,  12]))


## TSNE进行数据的降维,降维到2维空间中
tsne = TSNE(n_components = 2,perplexity =100,
            early_exaggeration =5,random_state=123) 

## 获取降维后的数据
tsne_outlier = tsne.fit_transform(Xs)
print(tsne_outlier.shape)
## 可视化在二维空间中的分布
plt.figure(figsize=(10,7))
sns.scatterplot(x=tsne_outlier[:,0], y=tsne_outlier[:,1],
                style = Y,hue = Y,s = 80)
plt.legend(loc = 1)
plt.grid()
plt.title("TSNE降维后数据的分布")
plt.show()

(672, 2)


## 介绍使用SVM分类器建立分类模型（这里只介绍这一种方法，其它的就不再介绍了）
## 建立非线性SVM模型，使用rbf核
rbfsvm = SVC(kernel  = "rbf",gamma=0.1, ## RBF核和对应的参数
             random_state= 1,C=10)
## 训练模型
rbfsvm.fit(X_train,y_train)
## 计算在训练集和测试集上的预测精度
rbfsvm_lab = rbfsvm.predict(X_train)
rbfsvm_pre = rbfsvm.predict(X_test)
print("训练集预测精度:",accuracy_score(y_train,rbfsvm_lab))
print("测试集预测精度:",accuracy_score(y_test,rbfsvm_pre))
print("训练集预测精度f1_score:",f1_score(y_train,rbfsvm_lab))
print("测试集预测精度f1_score:",f1_score(y_test,rbfsvm_pre))
print("在训练集上是否为异常值判断正确的精度为:\n",classification_report(y_train,rbfsvm_lab))
print("在测试集上是否为异常值判断正确的精度为:\n",classification_report(y_test,rbfsvm_pre))

## 通过参数调整，可能也可以获得较好的预测效果（也可能不会）

训练集预测精度: 0.996031746031746
测试集预测精度: 0.9285714285714286
训练集预测精度f1_score: 0.9696969696969697
测试集预测精度f1_score: 0.45454545454545453
在训练集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       1.00      1.00      1.00       471
           1       0.97      0.97      0.97        33

    accuracy                           1.00       504
   macro avg       0.98      0.98      0.98       504
weighted avg       1.00      1.00      1.00       504

在测试集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.96      0.97      0.96       156
           1       0.50      0.42      0.45        12

    accuracy                           0.93       168
   macro avg       0.73      0.69      0.71       168
weighted avg       0.92      0.93      0.93       168


## 对数据特征进行一次降维式的特征提取，看能否获得较好的结果。以PCA为例
from sklearn.decomposition import PCA

## 使用主成分分析对酒数据集进行降维
pca = PCA(n_components = 30,random_state = 123)
pca.fit(Xs)
## 可视化主成分分析的解释方差得分
exvar = pca.explained_variance_
plt.figure(figsize=(10,6))
plt.plot(exvar,"r-o")
plt.xlabel("特征数量")
plt.ylabel("解释方差大小")
plt.title("主成分分析")
plt.show()

## 先使用前10个主成分


pca_Xs = pca.transform(Xs)[:,0:10]
print(pca_Xs.shape)

## 数据切分训练集和验证机
X_train, X_test, y_train, y_test = train_test_split(pca_Xs, Y,test_size=0.25,random_state=2)
print("训练数据:",X_train.shape)
print("训练数据中异常值数量:",np.unique(y_train,return_counts = True))
print("测试数据:",X_test.shape)
print("测试数据中异常值数量:",np.unique(y_test,return_counts = True))

(672, 10)
训练数据: (504, 10)
训练数据中异常值数量: (array([0, 1], dtype=int32), array([471,  33]))
测试数据: (168, 10)
测试数据中异常值数量: (array([0, 1], dtype=int32), array([156,  12]))


## 获取降维后的数据可视化
## 可视化在二维空间中的分布
plt.figure(figsize=(10,7))
sns.scatterplot(x=pca_Xs[:,0], y=pca_Xs[:,1],
                style = Y,hue = Y,s = 80)
plt.legend(loc = 1)
plt.grid()
plt.title("PCA降维后数据的分布")
plt.show()


# 建立非线性SVM模型，使用rbf核
rbfsvm = SVC(kernel  = "rbf",gamma=0.1, ## RBF核和对应的参数
             random_state= 1,C=100)
## 训练模型
rbfsvm.fit(X_train,y_train)
## 计算在训练集和测试集上的预测精度
rbfsvm_lab = rbfsvm.predict(X_train)
rbfsvm_pre = rbfsvm.predict(X_test)
print("训练集预测精度:",accuracy_score(y_train,rbfsvm_lab))
print("测试集预测精度:",accuracy_score(y_test,rbfsvm_pre))
print("训练集预测精度f1_score:",f1_score(y_train,rbfsvm_lab))
print("测试集预测精度f1_score:",f1_score(y_test,rbfsvm_pre))
print("在训练集上是否为异常值判断正确的精度为:\n",classification_report(y_train,rbfsvm_lab))
print("在测试集上是否为异常值判断正确的精度为:\n",classification_report(y_test,rbfsvm_pre))

## 通过参数调整，可能也可以获得较好的预测效果（也可能不会）(利用主成分特征，效果有所增加)

训练集预测精度: 1.0
测试集预测精度: 0.9523809523809523
训练集预测精度f1_score: 1.0
测试集预测精度f1_score: 0.6666666666666666
在训练集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       1.00      1.00      1.00       471
           1       1.00      1.00      1.00        33

    accuracy                           1.00       504
   macro avg       1.00      1.00      1.00       504
weighted avg       1.00      1.00      1.00       504

在测试集上是否为异常值判断正确的精度为:
               precision    recall  f1-score   support

           0       0.97      0.97      0.97       156
           1       0.67      0.67      0.67        12

    accuracy                           0.95       168
   macro avg       0.82      0.82      0.82       168
weighted avg       0.95      0.95      0.95       168


## 使用参数搜索的基本模式
## 定义模型
svc = SVC(kernel  = "rbf")

## 定义网格搜索的参数
gammas=[0.001,0.001,0.1,10]
Cs = [0.5,5,50,100,1000]
para_grid = {"gamma": gammas,"C":Cs}

## 使用5折交叉验证进行搜索,使用均方根误差的负数作为得分
gs_svr = GridSearchCV(svc,para_grid,cv = 4,n_jobs=4,
                      scoring = "f1")
gs_svr.fit(X_train,y_train)

GridSearchCV(cv=4, estimator=SVC(), n_jobs=4,
             param_grid={'C': [0.5, 5, 50, 100, 1000],
                         'gamma': [0.001, 0.001, 0.1, 10]},
             scoring='f1')


print("最好模型使用的参数为:\n",gs_svr.best_params_)
rbfsvr_lab = gs_svr.best_estimator_.predict(X_train)
rbfsvr_pre = gs_svr.best_estimator_.predict(X_test)
print("训练数据集上的均方根误差:",f1_score(y_train,rbfsvr_lab))
print("测试数据集上的均方根误差:",f1_score(y_test,rbfsvr_pre))

## 可以发现在测试集上的预测精度更高了

最好模型使用的参数为:
 {'C': 5, 'gamma': 0.1}
训练数据集上的均方根误差: 0.7272727272727273
测试数据集上的均方根误差: 0.47058823529411764


## 将上面的程序进行相应的修改、组合、拼接，可以同时兼顾问题2中需要考虑的两点。

	时间	基站编号	小区编号	本地小区标识	上行可用的PRB个数	下行可用的PRB个数	上行PhysicalResourceBlock被使用的平均个数	下行PhysicalResourceBlock被使用的平均个数	上行PUSCH的PhysicalResourceBlock被使用的平均个数	小区内的平均用户数	...	MR测量上报RSRP在Index0区间的次数	MR测量上报RSRP在Index1区间的次数	MR测量上报RSRP在Index2区间的次数	MR测量上报RSRP在Index3区间的次数	MR测量上报RSRP在Index4区间的次数	小区QCI为1的DRB业务PDCPSDU上行期望收到的总包数	小区QCI为1的DRB业务PDCPSDU下行空口发送的总包数
0	2021-08-28 00:00	1200071	26019014	0	100	100	6.4553	31.5065	2.1503	22.6341	...	3	16	17	178	9943	2218	2179
1	2021-08-28 00:00	1200071	26019015	1	100	100	8.8281	21.5452	4.2220	19.1064	...	2	4	21	180	8427	1510	1423
2	2021-08-28 00:00	1200071	26019016	2	100	100	5.5024	11.9276	1.8618	15.4729	...	2	4	48	575	6371	2030	1919
3	2021-08-28 00:00	1200071	26019023	3	100	100	3.1113	0.9637	0.0077	0.0509	...	0	0	0	3	31	0	0
4	2021-08-28 00:00	1200071	26019024	4	100	100	2.6004	1.0051	0.1201	0.6355	...	0	1	1	9	285	7526	7780

	时间	基站编号	小区编号	本地小区标识	number
0	2021-08-28	1200071	26019014	0	22.6341
1	2021-08-28	1200071	26019015	1	19.1064
2	2021-08-28	1200071	26019016	2	15.4729
3	2021-08-28	1200071	26019023	3	0.0509
4	2021-08-28	1200071	26019024	4	0.6355

	时间	基站编号	小区编号	本地小区标识	number	outlier
0	2021-08-28	1200071	26019014	0	22.6341	None
1	2021-08-28	1200071	26019015	1	19.1064	None
2	2021-08-28	1200071	26019016	2	15.4729	None
3	2021-08-28	1200071	26019023	3	0.0509	None
4	2021-08-28	1200071	26019024	4	0.6355	None

	时间	基站编号	小区编号	本地小区标识	number	outlier
0	2021-08-28 00:00:00	1200071	26019014	0	22.6341	False
1	2021-08-28 00:00:00	1200071	26019015	1	19.1064	False
2	2021-08-28 00:00:00	1200071	26019016	2	15.4729	False
3	2021-08-28 00:00:00	1200071	26019023	3	0.0509	False
4	2021-08-28 00:00:00	1200071	26019024	4	0.6355	False
...	...	...	...	...	...	...
40363	2021-09-25 23:00:00	1200075	26019033	7	7.0792	False
40364	2021-09-25 23:00:00	1200075	26019034	8	1.3647	False
40365	2021-09-25 23:00:00	1200075	26019035	9	2.4111	False
40366	2021-09-25 23:00:00	1200075	26019027	10	11.7958	False
40367	2021-09-25 23:00:00	1200075	26019028	11	15.4725	False

	时间	number	outlier
18	2021-08-28 00:00:00	15.5100	0
76	2021-08-28 01:00:00	13.5003	0
134	2021-08-28 02:00:00	9.3439	0
192	2021-08-28 03:00:00	9.3711	0
250	2021-08-28 04:00:00	8.6753	0
...	...	...	...
40096	2021-09-25 19:00:00	29.3806	0
40154	2021-09-25 20:00:00	21.1153	0
40212	2021-09-25 21:00:00	20.7158	0
40270	2021-09-25 22:00:00	19.1725	0
40328	2021-09-25 23:00:00	16.9786	0

为所有的小区建立模型，输出其是否为异常值，使用的数据类似usenumdf.head()¶

利用modeldf的预测结果，建立一个预测未来是否为异常值的模型，二分类模型，由于小区之间的差异性较大，感觉单独为每个小区建立一个模型会更好一些。¶

XGBOD 效果较好，因为它使用了类别标签¶

二分类算法的效果会更好（针对类别不平衡的数据，还可以考虑使用一些数据平衡的方法）¶

针对某小区使用更多的特征¶

新的数据准备好后，可以采用与上面相同的方式，建立分类模型¶

注意：因为我们的小区较多，所以不可能手动的进行对每个小区，调出一组合适的建模参数，因此可以采用参数搜索的方式，自动的为每个小区建模获得较优的一组参数和模型¶

	上行可用的PRB个数	下行可用的PRB个数	上行PhysicalResourceBlock被使用的平均个数	下行PhysicalResourceBlock被使用的平均个数	上行PUSCH的PhysicalResourceBlock被使用的平均个数	ERAB建立成功总次数	ERAB建立尝试总次数	ERAB异常释放总次数	ERAB正常释放总次数	系统间切换出ERAB正常释放总次数	eNodeB内同频切换出成功次数	eNodeB间同频切换出成功次数	eNodeB内同频切换出执行次数	eNodeB间同频切换出执行次数	eNodeB内异频切换出成功次数	eNodeB间异频切换出成功次数
18	100	100	8.6102	21.4134	4.5893	3961	3961	0	4008	0	119	321	119	322	4	7
76	100	100	7.3371	17.7044	3.3577	4345	4345	0	4362	0	95	133	95	133	0	3
134	100	100	6.4366	7.5928	2.4736	3508	3508	0	3541	0	104	99	104	99	0	2
192	100	100	6.1496	7.7304	2.1808	3714	3716	0	3723	0	49	84	49	84	0	3
250	100	100	5.3418	5.1037	1.4129	3754	3754	0	3751	0	43	64	43	65	0	2
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
40096	100	100	16.6807	34.8728	12.2390	4821	4823	2	4582	0	868	1432	868	1433	1	71
40154	100	100	16.9814	25.6097	12.6263	3158	3161	0	3034	0	664	1396	667	1396	2	68
40212	100	100	13.6847	35.7595	9.5000	3082	3082	0	2992	0	561	915	561	915	0	45
40270	100	100	16.0122	37.4606	11.9351	2601	2603	0	2576	0	474	774	474	774	0	29
40328	100	100	8.1208	15.2549	4.3681	2553	2553	1	2525	0	308	470	308	470	0	21